Graphics Plus/Graphics Plus.iso/general/procssng/ccs/ccs-11tl.lha/lbl/hips/sources/vfft/realvfft.cxx

Text File | 1992-03-12 | 4.8 KB | 169 lines

/* realvfft.cxx
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

This software is copyright (C) by the Lawrence Berkeley Laboratory.
Permission is granted to reproduce this software for non-commercial
purposes provided that this notice is left intact.

It is acknowledged that the U.S. Government has rights to this software
under Contract DE-AC03-765F00098 between the U.S.  Department of Energy
and the University of California.

This software is provided as a professional and academic contribution
for joint exchange. Thus, it is experimental, and is provided ``as is'',
with no warranties of any kind whatsoever, no support, no promise of
updates, or printed documentation. By using this software, you
acknowledge that the Lawrence Berkeley Laboratory and Regents of the
University of California shall have no liability with respect to the
infringement of other copyrights by any part of this software.

For further information about this notice, contact William Johnston,
Bld. 50B, Rm. 2239, Lawrence Berkeley Laboratory, Berkeley, CA, 94720.
([e-mail])

For further information about this software, contact:
    Jin Guojun
    Bld. 50B, Rm. 2275, Lawrence Berkeley Laboratory, Berkeley, CA, 94720.
    [e-mail]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/

dimen2size = dimen1len * row << 1; /* use non complex represent complex size */

if (dflag)    {
double    *itemp,
    *ibuf = nzalloc(fsize, sizeof(*ibuf), "fibuf"),
    *obuf = nzalloc(dimen2size*(dimens?1:frm), sizeof(*obuf), "fobuf");

    uimg.o_form = dimens ? IFMT_DVFFT2D : IFMT_DVFFT3D;
    uimg.pxl_out = 16;
    (*uimg.header_handle)(HEADER_WRITE, &uimg, argc, argv, True);

    if (uimg.pxl_in < sizeof(*itemp))
        itemp = nzalloc(fsize, uimg.pxl_in, "itemp");
    else    itemp = ibuf;

    DBwork_space_init(MAX(MAX(cln, row), frm));

    for (f=0; f<frm; f++){
        i = upread(itemp, uimg.pxl_in, fsize, stdin);
        if (i != fsize)
        syserr("frame %d [%ld] read %d", f, fsize, i);
        if (uimg.in_form != IFMT_DOUBLE){
        switch(uimg.in_form){
        case IFMT_BYTE:    btod(itemp, ibuf, fsize);
            break;
        case IFMT_SHORT:    stod(itemp, ibuf, fsize);
            break;
        case IFMT_LONG:    itod(itemp, ibuf, fsize);
            break;
        case IFMT_FLOAT:    ftod(itemp, ibuf, fsize);
        }
        }

        DBvfft2d(ibuf, cln, row, obuf + (dimens ? 0 : f*dimen2size));

        if (dimens){
        i = fwrite(obuf, sizeof(*obuf), dimen2size, out_fp);
        if (i != dimen2size)
            syserr("2d[%d] write %d", dimen2size, i);
        }
    }
    if (!dimens && frm>1){
    int    r;
    DBCOMPLEX    *cp;
    register DBCOMPLEX    *cmptr;
    register int    c, n=frm, dimen1size = dimen1len;

        load_DBw(frm);

        for (r=0; r<row; r++)
            for (c=0; c<dimen1size; c++){
            cp = (DBCOMPLEX*)obuf + r*dimen1size + c;
            cmptr = DBcmpx_in;
            for (i = 0; i < n; i++, cp+=dimen2size>>1)
                cmptr[i] = *cp;

            DBFourier(cmptr, n, DBcmpx_out);

            cp = (DBCOMPLEX*)obuf + r*dimen1size + c;/* store back to dst */
            cmptr = DBcmpx_out;
            for (i = 0; i < n; i++, cp+=dimen2size>>1)
                *cp = cmptr[i];
            }
    /*    row = dimen1size;    real columns for IFMT_COMPLEX */
        fsize = frm*dimen2size;
        i = fwrite(obuf, sizeof(*obuf), fsize, out_fp);
        if (i != fsize)
            syserr("3d[%d] write %d", fsize, i);
    }
}
else    {
float    *itemp,
    *ibuf = nzalloc(fsize, sizeof(*ibuf), "fibuf"),
    *obuf = nzalloc(dimen2size*(dimens?1:frm), sizeof(*obuf), "fobuf");

    uimg.o_form = dimens ? IFMT_VFFT2D : IFMT_VFFT3D;
    uimg.pxl_out = 8;
    (*uimg.header_handle)(HEADER_WRITE, &uimg, argc, argv, True);

    if (uimg.pxl_in < sizeof(*itemp))
        itemp = nzalloc(fsize, uimg.pxl_in, "itemp");
    else    itemp = ibuf;

    work_space_init(MAX(MAX(cln, row), frm));

    for (f=0; f<frm; f++)    {
        i = upread(itemp, uimg.pxl_in, fsize, stdin);
        if (i != fsize)
        syserr("frame %d [%ld] read %d", f, fsize, i);
        if (uimg.in_form != IFMT_FLOAT){
        switch (uimg.in_form)    {
        case IFMT_BYTE:    btof(itemp, ibuf, fsize);
            break;
        case IFMT_SHORT:    stof(itemp, ibuf, fsize);
            break;
        case IFMT_LONG:    itof(itemp, ibuf, fsize);
            break;
        }
        }

        vfft2d(ibuf, cln, row, obuf + (dimens ? 0 : f*dimen2size));

        if (dimens){
        i = fwrite(obuf, sizeof(*obuf), dimen2size, out_fp);
        if (i != dimen2size)
            syserr("2d[%d] write %d", dimen2size, i);
        }
    }
    if (!dimens && frm>1){
    int    r;
    COMPLEX    *cp;
    register COMPLEX    *cmptr;
    register int    c, n=frm, dimen1size = dimen1len;

        load_w(frm);

        for (r=0; r<row; r++)
            for (c=0; c<dimen1size; c++){
            cp = (COMPLEX*)obuf + r*dimen1size + c;
            cmptr = cmpx_in;
            for (i = 0; i < n; i++, cp+=dimen2size>>1)
                cmptr[i] = *cp;

            Fourier(cmptr, n, cmpx_out);

            cp = (COMPLEX*)obuf + r*dimen1size + c;/* store back to dst */
            cmptr = cmpx_out;
            for (i = 0; i < n; i++, cp+=dimen2size>>1)
                *cp = cmptr[i];
            }
    /*    row = dimen1size;    real columns for IFMT_COMPLEX */
        fsize = frm*dimen2size;
        i = fwrite(obuf, sizeof(*obuf), fsize, out_fp);
        if (i != fsize)
            syserr("3d[%d] write %d", fsize, i);
    }
}